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We present a history-dependent Monte Carlo scheme for the efficient calculation of the free- 
energy of quantum systems, inspired by the Wang-Landau sampling and metadynamics method. 
When embedded in a path integral formulation, it is of general applicability to a large variety 
of Hamiltonians. In the two-dimensional quantum Ising model, chosen here for illustration, the 
accuracy of free energy, critical temperature, and specific heat is demonstrated as a function of 
simulation time, and successfully compared with the best available approaches, particularly the 
Wang-Landau method over two different Monte Carlo procedures. 
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Calculating certain thermodynamical quantities, such 
as the free-energy (FE) or the entropy, by Monte Carlo 
(MC) simulation is a notorious difficult problem. The 
difficulty arises because standard MC [l[ is devised so 
as to generate configurations A distributed according to 
their Boltzmann weight Px = e~ l3Ex /Z, where Ex is the 
energy of the configuration A and Z = J^x s~ /3Ex the 
partition function. This is efficient if we are interested 
in calculating quantities like the average energy (E) = 
Ylx ExPx, since the configurations generated by MC 
are just those that contribute significantly to the average. 
Calculating, however, the free-energy, F — —ft~ 1 log(Z), 
requires a knowledge of the partition function Z which is 
not accurately given by the simulation. 

A major step forward, in this respect, came with the 
Wang-Landau (WL) idea . In a nutshell, since: 



Z = e' pEx = J dEg(E)e 
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where g(E) — J2x^(^ ~ Ex) is the density of states 
with energy E, if we devise a MC that generates config- 
urations distributed according to Y/g(Ex), then we will 
effectively reconstruct the full histogram for g(E) in a 
single simulation. This allows computing the partition 
function Z, and hence all thermodynamical quantities, 
at any temperature T = l/feg/3, where ks is the Boltz- 
mann constant. This is particularly useful if the system 
can undergo a first-order phase transition. Indeed, using 
the WL approach, the system can diffuse over barriers 
between different local minima following pathways that 
would represent, in normal finite-T MC, "rare events". 

This discussion applies to classical systems; How 
should one proceed for a quantum system? Q Consider, 
to fix ideas, the transverse-field quantum Ising model 
(QIM): 
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were of, af are Pauli matrices, J > is an exchange con- 
stant, h and T are respectively the longitudinal and trans- 
verse magnetic field, and (ij) denotes nearest-neighbors 
on a lattice of N sites. The partition sum Zqim = 
(A|e~^ H |A), where X — {<Tj=i...jv} is a configura- 
tion of all N spins, involves now a matrix element of 
e —pH rpjjg g rs ^. s |- e p towards rewriting it in a form sim- 
ilar to Eq. ([1]) consists in performing a Suzuki- Trotter 
decomposition [4| , leading to a path-integral expression 
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Effectively, we have a classical system with an extra 
time dimension, whose configurations A, over which we 
sum, are given by X = {<7j =1 ...jy.p =1 ...p}. The extra in- 
dex p labels the P Trotter slices in the time direction [f| . 
In the QIM case, the action A reads: 



A{X) = N 



JU T +J ± K Y - hM-x - -x In C 
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where U Y = ~(NP) 1 E p Efe-j <t <,p (T 3,p is thc 
classical interaction energy per spin, K^- = 
— (NP)~ 1 J2i P a i,p cr i,p+i is the quantum "kinetic 
energy" per spin, M^- = (JVP) -1 ^,,^ is the mag- 
netization per spin, J 1 - = -(P/2/?)ln [tanh (/3T/P)} > 
is the ferromagnetic coupling between adjacent spins 
in the time-direction, and C 2 = (1/2) sinh (2/3T/P). 
By introducing a multi-dimensional density of states 

g(U, K,M) = Ex S ( U - U x) S ( K - K x)k M ~ M x) we 
can easily rewrite: 

Zqim « JdUdKdM g(U,K,M) e -PA(u,K,M) 

= J dUdKdM e-Pnu,*,*) , ( 5 ) 

where A(U, K, M) = N [JU + J X K - hM - (P//3) In C] , 
and F(U, K, M) defines the FE as a function of 
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(U,K,M). For h — the relevant coordinates are 
two, U and K. Using the WL idea to reconstruct 
Zqim for all values of /3 and V requires now sampling a 
two-dimensional density of states histogram g(U,K) in 
terms of which Zqim oc JdUdKg(U, K )e-P N ^ U+J±K \ 
This approach is, however, not very efficient (see below). 

A much more convenient ("state-of-the-art") route is 
based on the so-called stochastic series expansion (SSE) 
0, 0] , and involves using a WL approach to reconstruct 
the coefficients g(n) — Tr(—H) n of a high-temperature 
expansion of the partition function Z = ^2 n ((3 n /n\)g(n) 
[§]. The SSE approach is particularly suited to treat 
quantum spin systems and other lattice quantum prob- 
lems, but is in general not straightforward, for instance, 
for quantum problems on the continuum. 

We propose here a new method to effectively calculate 
the FE of a quantum system. Our approach is based on a 
path-integral formulation and can be easily extended to 
complicated off-lattice quantum problems. The crucial 
ingredients were borrowed from the WL method and the 
metadynamics approach, a method which proved useful 
for exploring the FE landscape of complex classical sys- 
tems [9] as a function of many collective variables (CVs) 
S = (Si, • • • Sd)- 

In metadynamics, sampling is enhanced introducing a 
history-dependent potential V G (S,t), defined as a sum 
of Gaussians centered along the "walk" in CVs-space, 
that in time "flattens" the FE histogram as a function 
of the CVs: V G (S,t -v oo) F(S) [10]. This ap- 
proach has been mainly used within molecular dynam- 
ics. During the simulation the system is guided by the 
action of two forces, the thermodynamic one, which move 
it towards the local FE minimum, and that due to the 
history-dependent potential, which pushes it away from 
local minima. 

We show here how to integrate metadynamics in a MC 
procedure, in particular in a path-integral MC (PIMC), 
to sample the FE landscape of quantum systems as a 
function of physically relevant CVs. Again, we illustrate 
this approach in the quantum Ising model where we re- 
construct the FE as a function of three CVs, the mag- 
netization M, the potential energy U and the kinetic en- 
ergy K. As we will show, a calculation performed at a 
single point (J3,T,K) in parameter space, is sufficient to 
obtain the FE in a whole neighborhood of that point. 
The method is tested by comparing its efficiency against 
the state-of-the-art WL-SSE method [8j, or a WL over 
a standard PIMC Q: we prove that our approach is at 
least as good as the WL-SSE on a lattice quantum prob- 
lem, as well as being physically transparent and easily 
generalizable to different models. 

Given the classical-like path-integral expression for the 
partition function of our quantum model, for instance 
Z ~ Y/x e ~^ > see Eq. (|3J), we first define a small 
number d of CVs S/(A), I — l...d, which appear in 
the action A(X) = A(S(X)): in the QIM case there 
are d — 3 physically meaningful CVs, the potential en- 



ergy 5*1 = U, the kinetic energy Sa = K and the mag- 
netization S3 = M, in terms of which the action is 
A(S) = N[JU + J^K - hM - (P//3) In C] . Next, we 
perform a Metropolis walk in configuration space {A} 

in which the transition probability from A to A is mod- 
ified adding to the action a history-dependent potential 
V G (S(X),t): 
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where 5 A — A{X ) — A(X) is the change in action and 
5V G (t) = V G {S(x'),t) - V G (S(A),i). Whether or not 
a move is accepted, we update V G by adding to it a 
small localized repulsive potential (a Gaussian in nor- 
mal metadynamics 9]). Technically, this is best done 
by grid-discretizing the CVs-space and keeping track of 
V G (S( k \t) only at grid points S^; the value of V G 
at a generic point S(A) is then calculated by a lin- 
ear interpolation C from the neighboring grid-values: 
Vb(S(A),i) = £(Vb(SW,i)) where £(...) is the linear 
interpolation function, and S^ k \ k = I... 2d, are the 
points of the grid nearest-neighbors of S(A). In this 
scheme, the potential V G is updated on the neighboring 
grid-points S( fe ) as: 
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V G (S^,t+l) = V G (S^,t)+wJ[ I —— zc 1 

where the (+) sign is used if S^ < Si(X) and the (— ) 
sign otherwise, AS; is the spacing of the grid in the Si di- 
rection and w is a parameter that determines the speed of 
the FE reconstruction. Therefore, like in WL, the accep- 
tance changes every time a move is accepted or rejected, 
and the "walk" in configuration space is intrinsically non- 
Mar kovian (it depends on the history). At the beginning 
of the simulation the potential V G (S^ k \ t = 0), stored on 
the grid, is set to zero. Then, as the system moves in 
configuration space, V G is updated at each move as in 
Eq. ([7]). After a sufficient time, V G will approximately 
compensate the underlying FE profile flo| . A further im- 
provement can be obtained by taking as estimator of the 
FE not just a single profile V G , but the arithmetic aver- 
age of all the profiles between a "filling" time tp and the 
total simulation time t to t- 
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This reduces the error of the method, which drops fast 
to zero for large t to t — ^ [9j. 

When F(U, K, M) for a given value of the external pa- 
rameters (/?, T,h) is known, one can readily recalculate 
the new FE profile for a whole neighborhood in param- 
eter space. The equations for this extrapolation can be 
written as: 
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By logarithmic integration of F(U, K, M) with respect to 
one or more variables we immediately get the free-energy 
as a function of a reduced number of CVs. For instance: 



as with a direct application of WL to PIMC in which 
the two-dimensional g(U, K) is calculated. For the same 
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Fig. [Qshows F(M) for the QIM on a 8x8 lattice (N = 64 
spins), with P = 30 Trotter slices at two different points 
in parameter space. The agreement between the refer- 
ence F(M) and that calculated from F(U, K, M) is good, 
even if we extrapolate the F(U, K, M) from the ordered 
to the disordered side (or viceversa) of the phase transi- 
tion line. Thus with a single calculation of F(U, K, M) at 
a point (T, F, h) in parameter space, we can get reliable 
information for F(U, K, M) in a whole neighborhood of 
that point (see inset). 
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FIG. 1: Free energy profile for the 8x8 QIM, as a function 
of the magnetization for two different parameter values. The 
results at (fcsT = 1.86, T = 2.0, h = 0) (in units of J) are ob- 
tained by first calculating F(U, K, M), and then performing 
a logarithmic integration, Eq. (|12[) , to calculate F(M). The 
results at (fcgT = 2.2, T = 2.2, h — 0.02) are instead obtained 
by extrapolating the previous F(U, K, M) using Eqs. (I9I11|I . 
and then integrating to obtain F(M). As a reference for the 
comparison we use the results of an accurate umbrella sam- 
pling calculation [Tl[ (solid line). The inset shows the phase 
diagram of the model and the circle suggests the size of the 
extrapolation region. 

In order to test the efficiency of the proposed method 
we compare it with a SSE-WL simulation 0,0, as well 
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FIG. 2: Critical temperature for a 8 x 8 QIM, as a function of 
the MC time, calculated using three different methods: SSE- 
WL (solid circles) , our method (solid squared) and the PIMC- 
WL algorithm (inset). The reference (red line with error bars) 
is obtained by a long PIMC simulation. 

system of Fig. [T]we estimate T c (conventionally defined 
as the temperature at which the specific heat reaches its 
maximum value) as a function of the MC time with the 
three methods. The results are shown in Fig. [51 As a 
reference, we also computed T c by a very long PIMC cal- 
culation (red line with error bars in the Figure). In the 
SSE+WL calculation the histogram is considered "flat" 
when for all the values of n the histogram is larger than 
95 % of its average [l2j (the limit of 80 % suggested in 
Ref 0] leads, for this specific system, to systematic er- 
rors, data not shown). Instead, for PIMC-WL the 80 % 
limit is sufficient to reach convergence. The specific heat 
for our method was calculated computing a F(U,K) at 
k B T/J = 1.8 and T/J = 2.0, h/J = 0.0 and extrapo- 
lating in temperature according to Eq. ©. The grid 
spacing in the U and K directions was of 10 and 1 en- 
ergy levels respectively. However we needed a finer grid 
spacing of 1 also for U, for states with U < —1.79167, 
in order to avoid systematic errors that generally tend to 
arise close to the parameter boundary values. 

In order to extrapolate the FE in a meaningfull tem- 
perature interval AT ~ ± J/fc^ including the peak of the 
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specific heat, it is necessary to obtain quickly a large max- 
imum value of Vq ~ 80/cbT/ J for the system considered 
here. This is accomplished by starting the simulation 
with w = 8 • 1CT 3 decreasing it up to 1CT 4 in 2 • 10 6 MC 
steps (tp in Eq. ©), then w is not changed anymore, 
and the free energy is estimated using Eq. (jHJ). It is clear 
from the previous discussion that the optimal "filling" 
protocol is system-dependent. 

As shown in Fig. [2j using our approach we can obtain 
T c within the PIMC error bar, with an efficiency similar 
to the SSE+WL algorithm. The PIMC+WL method is, 
by comparison, an order of magnitude slower (Fig. [2j 
inset). Of course, the efficiency of the approach presented 
here is strongly influenced by the temperature where the 
reconstruction is performed, which should not be too far 
from T c (~ 10% smaller in the example considered here). 
However, T c can always be estimated approximately, e.g., 
by performing a preliminary calculation on a system of 
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FIG. 3: Specific heat as a function of the temperature for the 
32 x 32 QIM using two different methods, the PIMC technique 
(solid circles) and the proposed method (solid squares). The 
inset shows how the estimate of T c evolves as a function of 
the MC time using the present scheme. 

Fig. [3] shows the specific heat as a function of T for a 
larger system, N = 32 x 32, with P = 100 Trotter slices, 
calculated with PIMC and with the present method. Also 
in this case we computed F(U, K) and extrapolated in 
temperature according to Eq. with a grid spacing 
of 150 and 10 energy levels in U and K directions re- 
spectively ( no need to reduce the grid spacing near the 
parameter boundaries since the free energy is very high 
there). In this case w was decreased from 10 _1 to 5 ■ 10~ 3 
in 1.1 • 10 6 MC steps. After this time the free energy is 



estimated using Eq. {SJ. As shown in Fig. [3l our ap- 
proach reproduces the specific heat accurately between 1 
and 3k bT/ J. In the inset we show how T c converges as a 
function of the MC time. Remarkably, even for this much 
larger 32 x 32 system the convergence of T c needs roughly 
the same order of magnitude of MC steps of those needed 
for the small 8x8 system. 

In conclusion, we have introduced an efficient history- 
dependent Monte Carlo scheme that allows the accu- 
rate calculation of the free energy landscape of quantum 
systems. The proposed approach was tested on a two- 
dimensional quantum Ising model, where we reconstruct 
the free energy as a function of two and three collective 
variables. This allows reproducing the thermodynamic 
properties in a whole neighborhood of the point in pa- 
rameter space at which the calculation is performed. The 
number of MC steps that are necessary to estimate T c in a 
relatively large system (32 x 32 x 100) is of the same order 
as that required in a small system (8 x 8 x 30). The effi- 
ciency in estimating T c is similar to that of SSE+WL, the 
state-of-the-art approach. Based on path-integral MC, 
our method can however be directly applied to continu- 
ous, off-lattice quantum problems, where SSE would be 
harder to implement. 
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